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Three-dimensional simulations with fully resolved hydrodynamics are performed to study the col¬ 
lective motion of model swimmers in confinement. We show that certain swimming mechanisms 
can lead to traveling wave-like collective motion even without any direct alignment mechanism. It 
is also shown that by varying the swimming mechanism, this collective motion can be suppressed, 
contrary to the perception that hydrodynamic effects are completely screened at high volume frac¬ 
tion. From an analysis of bulk systems, it is shown that this traveling wave-like motion, which can 
be characterized as a pseudo-acoustic mode, is mainly due to the intrinsic swimming property of 
the particles. 
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Active matter systems such as groups of microorgan¬ 
isms, birds or fish often form swarms or flocks, and typ¬ 
ically show complex collective behavior BH3|. Of special 
interest for scientific study are active microscopic sys¬ 
tems, e.g. microorganisms or Janus particles, because 
they are particularly well suited for experiments. Under¬ 
standing the collective dynamic behavior of these systems 
can be useful for both science and industrial applications, 
for example, to explain the formation of biofilms [4j or to 
design targeted drug delivery systems [5J. Such systems 
typically exist under some type of confinement, which has 
been shown to exert a strong influence on the collective 
behavior hud]. So far, most analysis of such systems 
have either neglected hydrodynamic interactions or only 
considered them in a far field approximation HEnnioca, 
even though it is known that the hydrodynamic interac¬ 
tions can dramatically affect the physical properties of 
the systems m- Theoretically, it’s impossible to take 
into account the full hydrodynamic interactions in many 
body systems, and numerically it’s very computationally 
expensive. While several studies have succeeded in in¬ 
cluding hydrodynamic interactions HMD], none have 
focused on the collective dynamic behavior of 3D ac¬ 
tive systems under confinement. The recent work by 
Zottl and Stark[l7j is related to ours, but deals with 
pseudo-2D systems under strong confinement (small sep¬ 
aration of walls), high volume fraction and strong swim¬ 
ming strength, which are quite different conditions from 
those of the present work. Not surprisingly, we are able 
to observe quite different dynamic behavior. 

We have conducted direct numerical calculations for 
3D systems of spherical swimmers in a host viscous fluid 
with fully resolved hydrodynamics. Under the confine¬ 
ment of two flat parallel walls, we find collective motion 
of swimmers with traveling wave characteristics, which 
propagates back and forth between the walls. Although 
similar traveling wave-like motion has already been ob¬ 
served in systems with explicit alignment interactions be¬ 
tween the particles ii, here we show that hydrody¬ 
namic interactions alone are enough to exhibit this be¬ 


havior. We also confirmed that this traveling wave-like 
behavior is not due to the presence of the walls, but is 
just a manifestation of the bulk behavior. Crucial to ob¬ 
serve this behavior are the strength and type (pushers vs. 
pullers) of swimming, not the degree of global alignment. 

As the numerical model of microswimmers, the 
Squirmer Model was employed. In this model, a self- 
propelled particle is modeled as a spherical object with a 
modified stick boundary condition at its surface EDI El]: 

00 2 

u s {r) = ^2———B n P^(cos9)smdd, (1) 

z ' nln + 1 

2 = 1 X ' 

where f is a unit vector directed from the center of a 
squirmer to a point on its surface, 0 = cos -1 (r • e) the 
polar angle between r and the swimming direction e, 9 
the unit polar angle vector. P' n is the derivative of the 
Legendre polynomial of n-tlr order, and B n is the magni¬ 
tude of each mode. Here, the radial deformation has been 
ignored, so this surface velocity has only the tangential 
component, and is responsible for the self-propulsion of 
the swimmers. 

In this work, only the first two modes in Eq. 0 were 
retained: 

u s {9) = B\ ^sind + — sin 26^ 9. (2) 

The coefficient of the first mode, B\, determines the 
swimming velocity of an isolated squirmer, as Uq = 
2/3Bi. The coefficient of the second mode, B 2l deter¬ 
mines the strength of stirring, or the stresslet E2J • Using 
only the first two modes, the type of swimming can be 
modified by changing the sign of the second coefficient, a 
in Eq. 0. Negative values of a describe pushers, which 
generate an extensile flow field in the propeling direc¬ 
tion, and positive values describe pullers, which generate 
a contractile flow (Fig. [l]). Because these two swimming 
mechanisms result in distinct flow fields, different swim¬ 
mers can exhibit vastly different collective motion, as we 
will show in this paper. 
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FIG. 1: Schematic representation of pushers and pullers 
within the squirmer model. The swimming direction, and the 
fluid and surface flows generated by squirmers are shown. 

In this study, the squirmer model is incorporated into 
the Smoothed Profile Method (SPM) IMl!25| . which is 
an efficient calculation scheme for solid/fluid two phase 
dynamics problems with full hydrodynamics. In this 
method, the sharp boundary between the solid and fluid 
domains is not considered explicitly. Instead, a diffuse 
interface of finite width £ is introduced, and the solid 
domain is represented via an order parameter <j>, which 
takes a value of 1 in the solid domain, and 0 in the fluid 
domain. Using this continuous order parameter, all the 
physical quantities can be expressed as field variables on 
Cartesian grids, and the calculation cost can be reduced 
drastically. As the governing equations for the total fluid, 
we employ a modified incompressible Navier-Stokes equa¬ 
tion: 

Pi {dt + u ■ V) u = V ■ cr + Pi (0/ p + / sq ) , 

V • u = 0, (3) 

u = (1 — (f))U{ + 0Up, 

(/)Up = 'y ^ [V i -b 12^ x Ri \, 
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where u is the total velocity field, (1 — 4>) u i and 4 >u p 
the velocity field of fluid and particles, respectively. The 
position, velocity and angular velocity of particle i are 
expressed as R,, Vi and IV The stress tensor and the 
mass density of the host fluid are cr, p f . The body force 
generated from the rigidity condition of particles is calcu¬ 
lated as 4>f p , and / stands for the force field resulting 
from the squirming motion [26]. The time evolution of 
the particles follows from the Newton-Euler equations: 

Ri = Vi , Qi = skew(fU) ■ Q L , (4) 

M p Vi =Ff + Ff, I p tli = Nf , 

where denotes the rotation matrix, M p the mass, 
skew(fij) the skew symmetric matrix of fl,, and F^ and 
IVf are the hydrodynamic force and torque acting on 
particle i. The force due to direct interactions between 
particles, Ff, is taken to be a truncated Lennard-Jones 
potential, in order to ensure that the excluded volume 


constraint is satisfied. See refs [22H2E3 for more detailed 
information of this method. 

Using this computational model, we conducted simu¬ 
lations of squirmer dispersion systems in a three dimen¬ 
sional domain. Throughout this paper, the diameter a 
and the boundary thickness £ of the squirmers are 4A, 
2A, where A is the grid spacing and the unit length, and 
the parameter B\ in Eq. ([2]) is 0.375, respectively. The 
shear viscosity ?y and the fluid mass density pf are set 
to one, and the unit time is then expressed as pA 2 / p. 
The particle Reynolds number of an isolated squirmer 
Re 0 = piU^a/r] is set equal to one. The Reynolds num¬ 
ber in the dispersion will be different from Re 0 due to a 
decrease in the propelling velocity with increasing volume 
fraction. We have ignored any effects due to thermal fluc¬ 
tuations. No external force is exerted, and the systems 
are buoyancy free. 

We first calculated the dynamics of N = 2607 swim¬ 
mers in a 64A x 160A x 64A rectangular system (volume 
fraction of 13%). The computational domain is periodic 
in the x and z direction, but confined between two flat 
parallel walls along the y direction; at y = 1A and 159A. 
These walls are made of 2048 spherical particles which 
are pinned and unable to move or rotate from the initial 
positions and orientations. The wall particles are regu¬ 
larly spaced and are defined to have the same diameter 
as the squirmers. The initial configuration of swimmers 
was randomly determined with respect to both position 
and orientation. We conducted simulations for various 
a values in order to investigate the effect of a on the 
dynamics. We have calculated the time evolution of the 
local density p , polar order Q\ and velocity v\, as a func¬ 
tion of distance from the walls. Here, Q\ = (e 3/ ) p i ane and 
v\ = (Uy/Uo)plane, where e y and V y are the components 
of e and V perpendicular to the walls, and (-(plane de¬ 
notes an ensemble average over all particles in parallel 
slabs of width 2A. Fig. [2] shows the results of the sys¬ 
tems with a = ±0.5 (see supplementary materials SI 
and S2). For pullers (a = 0.5), the time evolution of 
the density, polar order and velocity are shown, while for 
pushers (a = —0.5), only density is shown. Only the 
pullers show large density fluctuations, which we note 
are due to hydrodynamic interactions alone, as no direct 
alignment interaction has been included; pushers show 
no such behavior. First, the pullers build up two distinct 
traveling waves which bounce back and forth between 
the walls. Gradually, as one wave begins to dominate, 
they merge into a single propagating wave. We also con¬ 
sidered different wall separations, and found that such 
propagating wave-like motion can only be observed if the 
wall separation is larger than the characteristic size of 
this propagating flock. In our simulations, the threshold 
is W/a — 40, where W is the wall separation. In fig. [3j 
the traveling wave velocity is plotted for several different 
wall separations (the bulk sound velocity which will be 
explained later is also shown). They are quantitatively 
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(a) snapshots 



(b) local order (pullers) 



(c) local velocity (pullers) 



(d) density (pullers) 



(e) density (pushers) 
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FIG. 2: Simulation results for pullers and pushers (a = ±0.5) 
at volume fraction of 13% and wall separation W/a = 40. 
(a) Simulation snapshots at T = 11Tb, (b) local order of 
pullers defined as Q\ = (e v ) plane, (c) local velocity of pullers 
in y direction (perpendicular to the walls) defined as hi = 
(V y /Uo )plane (d,e) plane averaged density for pushers and 
pullers at each height y/W. The density is normalized by the 
global average value, pa, and time and velocity are normalized 
by To = W/Uo and Uq. where W represents the separation 
of parallel flat walls. The snapshots are of around t = UTo, 
which is marked by dashed line. 


the same and there is no clear separation dependence. 
From these measurements, this traveling wave-like mo¬ 
tion can be understood as an intrinsic property of the 
self-propelled systems, which is independent of the wall 
effects. 

To confirm this conjecture, further simulations were 
conducted in cubic systems with no confinement, un¬ 
der full periodic boundary conditions in all directions. 
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FIG. 3: The traveling wave velocity of confined systems as 
a function of wall separation W. The sound velocity of the 
bulk system is also shown. Here, only puller results are shown, 
as pushers don’t show wave-like behavior. All the values are 
normalized by Uq. 
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The linear dimension of the computational domain is 
L = 64A and the number of swimmers is 1043 ( volume 
fraction of 13 %, the same value as the previous simula¬ 
tions in the systems with walls). In order to study the 
density correlations, we calculated the dynamic structure 
factor S(k 1 Lo), which is just the Fourier transform (in 
both space and time) of the density-density correlation 
functionP?]: 


N 

Pk ( t ) = ^2 exp \-ik ■ R t (t)] (5) 

i =1 

1 f°° 1 

S(k,w) = — J — (p k {t)p- k (0))exp (iwf) dt (6) 

where k is the wave vector and w the angular frequency. 

The spherically averaged (over shells of wave vectors) 
structure factors are shown in Fig. [4] Shifted peaks were 
observed for both systems, pullers and pushers. These 
peak positions w s hift follow a Brillouin-like dispersion re¬ 
lation: 


k’shift — c| ^ | 


(7) 


where c is a coefficient which has the dimension of veloc¬ 
ity (in the case of normal Brillouin peaks, c represents the 
sound velocity). We can therefore say that these systems 
possess pseudo-acoustic mode dynamics, or a sound-like 
propagation mode for the density. No Rayleigh-like peaks 
were observed for small wave vectors, probably because 
the Brillouin-like peaks are much stronger (for wave vec¬ 
tors of large magnitude, Rayleigh-like peaks are recov¬ 
ered). The sound velocity of this pseudo-acoustic mode 
can be calculated from the dispersion relation,shown at 
the bottom of Fig. [4] We don’t see any difference between 
the sound velocities of pushers and pullers, although the 
peak intensities show very big differences. The calculated 
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value of the sound velocity for the pullers is quantita¬ 
tively the same as the propagating wave velocities in the 
confined systems, as shown in fig. [3j although the veloci¬ 
ties are slightly reduced for small wall separations (Sup¬ 
plementary Material S3). Therefore, we can conclude 
that the traveling wave-like behavior in the confined sys¬ 
tem is mainly due to the intrinsic pseudo-acoustic prop¬ 
erty of the swimmers themselves, and the wall effects are 
not crucial. 

Now we would like to discuss the origin of the pseudo¬ 
acoustic mode. From fig. [2j we can see that the particles 
show local polar order and move as a group collectively, 
even though individual particles move independently of 
each other (for individual particles motion, see Supple¬ 
mental Material S4). In puller systems, this cooperative 
swimming as a group can be understood as the origin of 
the pseudo-acoustic mode, or density propagation. Now 
let us consider a simplified set-up in which swimmers are 
placed along a straight line with a non-regular distribu¬ 
tion, and oriented parallel to this line, in such a way they 
will start swimming in the same direction (Supplemental 
Material S5 and S6). In this situation, due to the con¬ 
tractile flow along the swimming direction, pullers at¬ 
tract each other and tend to enhance the longitudinal 
mode of density fluctuations, resulting in the formation 
of flocks. On the other hand, pullers show the opposite 
tendency, which suppresses the longitudinal fluctuations. 
In the systems of interest in this paper, the situation 
is much more complicated. However, we can still say 
that the appearance of the large density fluctuations we 
present depends on the local interaction rather than on 
the global order. Actually, we never observed large fluc¬ 
tuations for pushers, even though the global polar order 
(Q = (lylEi^D) ii could be similar to pullers 
with a = 0.5. On the other hand, pullers with various 
values of a (0.3 < a < 0.9) do show traveling-wave like 
motion in confinement. In order to rule out the possibil¬ 
ity that this fluctuation is due to inertial effects, we also 
conducted a simulation with smaller value of Reynolds 
number (Re o = 0.1). The system showed the same col¬ 
lective behavior. We can therefore conclude that this 
behavior is purely hydrodynamic in origin. 

Several articles have shown similar results; some are 
consistent with our findings, while others are not. Simha 
and Ramaswamy predicted, using a continuum model 
that active matter suspensions with polar order can show 
sound like propagation modes [Til . They mainly studied 
the propagation mode of a director field, but also showed 
that the concentration field can exhibit the same type of 
propagation. Ezhilan et al. studied a similar continuum 
model [T2J. They reported pushers show stronger insta¬ 
bility in their concentration field, contrary to our results 
that pullers have much higher Brillouin-like peaks and 
are more likely to show density fluctuations. They em¬ 
ployed a far-held approximation, therefore they did not 
consider the finite size of the particles or the full hydro- 







FIG. 4: The calculation results for the spherically averaged 
Dynamic Structure Factor (DSF) S ( k,co ), where k = ^4 k is 
the wave vector, and uj the angular frequency (normalized by 
uq = 2n/to, where to = a/Uo). The top two figures show 
the simulation results for the DSF of (a) pullers at a = 0.5, 
and (b) pushers at a = —0.5. The bottom figure (c) shows 
the relation between wave vectors and the positions of the 
Brillouin peaks. Arrows in (a) and (b) show the positions of 
peaks. 


dynamic interactions. In our study, we took into account 
the full hydrodynamic interactions and finite volume of 
the swimmers. This difference is likely the cause for the 
contradictory results. Alarcon and Pagonabarraga[T5] 
reported the same kind of density fluctuations in bulk 
puller systems, using direct numerical simulations with 
the same model swimmer system. They focused on the 
flocking behavior in bulk systems, while we studied the 
behaviors in confinement, which makes the flocking more 
apparent, and also identified the pseudo-sound property 
of the systems. 
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In conclusion, we confirmed, using direct numerical 
simulations with full hydrodynamics, that swimmer dis¬ 
persions confined between flat parallel walls can exhibit a 
traveling wave-like motion. Furthermore, this non-trivial 
collective behavior is mainly due to the pseudo-acoustic 
properties of the swimmers themselves, with a purely 
hydrodynamic origin. This kind of phenomenon cannot 
be expected in conventional passive colloidal dispersions, 
since such density fluctuations would be suppressed due 
to viscous damping BUI . In future works, we will con¬ 
sider complex wall geometries, such as gear shaped ob¬ 
jects and circular systems jS]. Based on our current 
results, we expect that spherically shaped particles, with¬ 
out any apparent alignment interaction, can generate non 
trivial collective motion due solely to hydrodynamic in¬ 
teractions. 
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